%Written by Michael R. Springborn
%This version: Dec. 16, 2011

% Free to use and modify with proper acknowledgement. No warranties. 
%*******************************************************************************

%                         ***GENERAL DESCRIPTION***     

%The main purpose of this program is to set up and run MCMC simulations using the 
%Metropolis-Hastings algorithm to generate a sequence of random samples from the 
%(Bayesian) posterior distribution for the parameters of a logit or cauchit model. 

%The code in its present form is geared towards replicating the estimation 
%results in the paper 'Closing the Gap Between Risk Estimation and 
%Decision Making: Efficient Management of Trade Related Invasive Species Risk', 
%by Lieli and Springborn, forthcoming, Review of Economics and Statistics.

%Output is analyzed using code in "analysis.m" and results are reported in Tables
%3 and 4 of Lieli and Springborn. 

%Requires the custom function CosGibbsMH().  That function is parameterized with some 
%initial values using maximum likelihood estimation found in "runmle.m"

%*******************************************************************************
% 1. General
%*******************************************************************************

clear; clc;
%set the working directory.  
work_dir='D:\Bayesian_Analysis'
cd(work_dir)

%*******************************************************************************
% 2. The raw data 
%*******************************************************************************
importfile('D:\data\data_WRA.csv'); %For descriptions of field names see Excel file of same name in same folder.

%*******************************************************************************
% 3. Select options for the link model, simulation algorithm, and whether a subset of the data is to be used.
%*******************************************************************************

%Select the link model for modeling binomial outcomes: 
%model='logit'; %Logistic latent variable: p_i=logit(B0 + B1*X_i)
model='cauchit'; %Cauchy latent variable: p_i=logit(B0 + B1*X_i)

%Set parameters for the Metropolis-Hastings algorithm (Markov chain Monte Carlo procedure)
s=20000; %Iterations for the chain
M=11; %Number of chains

%The estimation may be run with all observations (N=370) or on a subsample (<370) for testing.  
% Generate 1,000 sets of random observation numbers for later use in randomly selecting subsets of the data set.
% Note: can run this random number generation once and then turn off to keep the same random selection.
    rand4sub=1; %set to 0 once a random observation numbers have been generated and saved.
    N=370; %370 observations in full data set
    if rand4sub
        ssmat=zeros(N,1000); for j=1:1000; ssmat(:,j)=randperm(370)'; end; 
        save D:\Bayesian_Analysis\subsampmat.mat ssmat
    end
    
%%If conducting estimation using the entire data set:
%1. comment out the for-loop over iternum below
%2. set iternum to empty:
iternum=[]; 

%%If conducting estimation using subsamples of the data set:
%1. remove commenting on the for-loop over iternum below
%2. set number of subsamples to evaluate:
IN=1000; 


%******************************************************************
%4. Run the MCMC simulation using the Metropolis-Hastings algorithm to generate a sequence of random samples from the posterior distribution.  
%******************************************************************
%Function -- CosGibbsMH()
    %implements the MCMC simulation... 
    %"Cos-": addresses the stratified sampling problem using the approach of Cosslett (1993)
    %"-Gibbs-": using Gibbs sampling of parameters (mu, sigma) for covariate X  [See Appendix B in Lieli and Springborn]
    %"-MH": using the Metropolis-Hastings algorithm

%Inputs to CosGibbsMH()
    %Qd: Assumed base rate of weeds (unconditional probability that Y=1).  Note: referred to in manuscript as "tau".
    %simp: Model for the density of the covariate sample, f(x;lambda) in Equation (13) of Lieli and Springborn.  [Also see Appendix B in Lieli and Springborn]
        %Reported results in Lieli and Springborn are for simp=1 (mean and variance of x set to their maximum likelihood estimates)
        %See "CosGibbsMH.m" for detail on all four variations of increasing complexity.
    %iternum: number for a particular random subset of the data.  Only used if running model on random subsets of the data.
    %m: chain index (will simulate several different Markov chains).  
    %s: iterations in each chain
    %Bin_DepVar,Score_assessor,Sc_Undes: data
    %model: link function (e.g. logit, cauchit).
    
for Qd=[0.02 .05] 
    for simp=1:4      
%        for iternum=1:IN
            for m=1:M
                CosGibbsMH(simp,Qd,m,s,Bin_DepVar,Score_assessor,Sc_Undes,iternum,model)
            end
%        end
    end
end


return
